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ABSTMCT 

The optimal control problem of a typical nuclear reactor 
power plant, which is described by a ninth-order nonlinear 
differential equation, having time-varying parameters, is 
considered. The nonlinear model complicates the optimal 
controller synthesis. Therefore, the approach of this \\fork 
is to approximate the response of the reactor system by 
that of a second-order linear model. The model parameters 
are chosen to minimize the derivations between the system 
and model responses using a search routine. The optimal 
feedback parameters computed for the second-order model is 
used for suboptimal control of the system. The model param- 
eters are updated to reflect the system nonlinearities as 
well as changes in the system parameters; the corresponding 
control scheme is adaptive. It is sliown that for the operat- 
ing conditions considered, the adaptive controller need not 
be on-line. 

Also, investigation of the effects of different weighting 
factors in the cost function, and the effect of various control 
rod configurations on the system response are presented. 
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INTRODUCTION 



I . 

Since the advent o£ nuclear reactor technology, safety 
is generally a foremost consideration in the design of a 
nuclear reactor and its control system. To attain the 
purpose of safe operation with some improvement in perfor- 
mance, there is a possible need to apply optimal control 
theory. A considerable amount of research effort has been 
directed toward the realization of optimal control systems, 
especially for such big-scale projects as the operation of 
a nuclear power reactor. It is reasonable to assume that 
nuclear pox'/er stations of the future x-/ill be operated 
entirely under the control of digital computers, which is 
accepted as the most suitable device to perform such a task. 

As mentioned by Sinha and Bereznai [Ref. 1], the over- 
simplified representation of the reactor dynamics, usually 
the one delayed neutron group model, is not satisfactory 
for control purposes, and their responses deviate consider- 
ably from that of the actual system. Therefore, the attempt 
is made to evaluate the response using a more complete 
representation of the reactor kinetics, and to include an 
adequate description of the controller mechanisms. 

In this paper, a realistic system of a nuclear reactor, 
xvhich is a ninth-order nonlinear and has time-varying 
parameters, is used. The measured quantity to be controlled 
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is the neutron or power level. Thougli any method that 
resolves the optimal control problem for this system is 
not known, it may be feasible to get a controller that 
will approximate the desired optimal response sufficiently 
close for practical purposes. A method that realizes 
such a suboptimal performance is described in this paper 
based on the work of Bereznai [Ref. 2] . It uses a second- 
order model to represent the reactor and the controller 
mechanism. The optimized parameters of the model are com- 
puted to approximate the response of the actual system 
using a search routine. Then, the feedback parameters 
required to use the suboptimal feedback controller for the 
nuclear reactor are easily computed. In this manner, the 
control system is adapted to com.pensate for the nonlinear 
system characteristics and for the changes in the system 
parameters . 
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SYSTEM REPRESENTATION 



II . 

T.he mathematical model for the plant under consideration 
is developed. To retain sufficient accuracy, the reactor 
kinetics represented by the six delayed neutron group is 
used. The unity feedback is used to study the closed- 
loop response of the system. 



A. REACTOR KINETICS 

The kinetic equations of a chain-reacting pile have 
been derived in the literature many times [Ref. 3] . The 
six delayed neutron group, space independent, source free 
model is described by the follovv’ing equations: 



dn (t) _ 6k(t) - B 
dt £ 



6 

+ Z X. C. (t) 
i = l " " 



( 2 . 1 ) 



dC.(t) 

dt 




n(t) 



X. 

1 



c. Ct) 



( 2 . 2 ) 



where i = 1, 2, ..., 6. 

All of these equations are linear in the so-called state 
variables, n and C^ (i = 1, ..., 6), and linear in tbe input, 
6k(t) , but the system is not jointly linear in state and 
control. Any feedback of n to formulate 6k(t) results in 
a nonlinear system. 
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B. TEMPERATURE REACTIVITY 



In most reactors, heat- transfer dynamics are coupled 
to neutron kinetics by a temperature coefficient of reactiv- 
ity which is generated by thermal expansion of the core and 
by change in various neutron cross sections. In other v;ords , 
thermal expansion reduces density, and thereby the modera- 
tion ability and reactivity. Also, the neutron capture-to- 
fission ratio of the fuel is altered by temperature changes, 
and a smaller effect, due to changes in absorption and 
fission cross sections, increases leakage and decreases 
reactivity. Further, thermal -reactivity contribution 
results from so-called Doppler broadening associated with 
the resonance region of the neutron energy spectrum. The 
temperature coefficient may usually be assumed to be a 
constant within the accuracy v/ith which it can be pre- 
dicted. Hence, the temperature reactivity may be expressed 
as : 

d6k. 6k a 

—rr^ = + — n(t) . (2.3) 

The significance of the temperature coefficient of reactivity 
is obvious from Figure 2.1. Even if a linear heat-exchange 
model is utilized, the temperature-reactivity feedback 
causes the reactor dynamical model to be nonlinear. The 
effective temperature is assumed to be proportional to the 
operating power level. 
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Figure (2.1) 

Feedback Loop of Temperature-Reactivity 
C. ABSORBER ROD 

The reactivity term 6k is the sum o£ the externally 
applied reactivity due to the absorber rod, and the change 
in the reactivity due to temperature variations. The mathe- 
matical representation o£ the reactivity due to the absorber 
rod is given in Re£erence 2: 



o 

II 

O 


/^W(V) 

0 


dt 




(2.4) 


Where W(V) = 0.02 


V(t) 


|V| 


< 15 


(2.5) 


= 0.3 




V 


> 15 




O 

II 




V 


< -15 





that is: 
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vct) 



Figure (2.2) 



Absorber Rod Characteristic Curve 



W(V) has the units mk/sec, and represents the reactivity- 
change due to the motion of the absorber rod. This formu- 
lation assumes that the control rod is near the center of 
its movement, vs^here the reactivity changes linearly \\'ith 
distance from the center of the core. 

D. DRIVE MOTOR 

The time response of the drive motor is given by: 



dV ^ 

3T T 




m 



T 



m 



( 2 . 6 ) 



v^?here G is gain, V is motor velocity 



e is power error and x is time 
^ m 

constant of motor. 
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E. OVERALL SYSTEM 



The block diagram of the nuclear reactor model considered 
in this thesis is represented in Figure (2.3). 




Figure (2.3) 

Block Diagram of Nuclear Reactor and Reactivity Mechanism 

In order to study the closed-loop response of the system, 
unity feedback is used. The error signal (e) that expresses 
the deviation between the demanded and actual power levels 
is amplified, and the output is applied to the absorber rod 
drive motor. The rate of reactivity insertion (6k^/sec) is 
proportional to the error signal until the motor reaches 
its maximum speed. The actual reactivity 6k, which is a 
measure of the criticality or multiplication factor of the 
reactor, is the difference between the reactivity worth of 
the absorber rod and the effect of the temperature change. 
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In state variable form the system equations that represent 
the reactor and the absorber rod mechanism in an open loop 



configuration are written in the following: 



fdn N 
dt 
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dt 
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dt 
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dt 
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0 
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dCg 
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-^6 
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d6k^ 




L. 
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-1 
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dt 
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-1 
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b 
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-i- 
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b 
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6k 
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( 2 . 7 ) 






/ 



Using matrix notation 



x=Ax+Be (2.8) 

and the system output is given by 

y = C^x (2.9) 

where = [1,0 ... 0] . 

The numerical values of the parameters are given in 
Appendix A. Hence, the parameters of reactor kinetics 
equations are typical values, and the parameters of the 
absorber rod mechanism are given in Reference 4. 
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III. DYNAMICAL BEHAVIOR OF SYSTEM 



The equations of the system described in Section II were 
solved on an IBM 360, Model 67 digital computer (W. R. Church 
Computer Center at the Naval Postgraduate School) using the 
method of digital simulation language (DSL) , which has been 
suitably modified to take into account the nonlinear nature 
of the differential equations. 

A. CLOSED-LOOP RESPONSE OF NUCLEAR REACTOR 

The closed-loop responses of the nuclear reactor to 
step-change in demanded power level, for various initial 
power level are shown in Figure (3.1) to Figure (3.8). 

The reference points such as the time to reach first full 
pov/er (1001 FP) , the amplitude and the time of maximum 
overshoot are tabulated in Table (3.1). The nonlinear 
characteristic of the system is clearly shown in the figures. 

B. EFFECT OF CONTROL-ROD TO THE RESPONSE OF NUCLEAR REACTOR 
Various control-rod configurations are tried to examine 

the effect to the system response. These are shown in 
Table (3.2). Hence, the number given to the control-rod 
shown in Table (3.2) are corresponding to the number of 
Figures (3.9) and (3.10). 

The configurations of 1, 2 and 3 have the same W 

m3- A. 

value and different slopes for the linear region. These 
system responses are shov.'n in Figure (3.9). In this case. 
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Comparison of Some Reference Points to Various Initial Power Levels 



Table (3.2) 

Control Rod Configurations 



No, 



Configuration 



Equation 






W 






W = O.OIV 


[vl < 30 


3 


0.3 






-30 








= 0.3 


V > 30 




l 

o 


30 

V 


o 

1 

li 


V < -30 




P 

II 

o 

o 

M 
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Figure (3.1) 

Steady State Response 




Figure (3.2) 

Step Response of Initial Power Level 951 
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Power Level (I) Power Level (%) 




Figure (3.3) 

Step Response o£ Initial Power Level 90 




Step Response of Initial Power Level 80"o 
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Figure (3.6) 

Step Response of Initial Power Level 50% 



Power Level Power Level 




Figure (3.7) 




Figure (3.8) 

Step Response of Initial Power Level 20 
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Figure (3.9) 

b G C 

Step Response of Reactor System Using Control Rods Designated 1, 2 and 3 






I 





100 




U 

O 

(/) 



LO 

C 



O 

ci 

C 

to 

• H 

LO 

O 

» 



o 

Pi 




K)^ 



to 

o ^ 

^.H 

-5 (O 

to>^ 

• H 

pH S 

o 

(/') 

>s 

CO 

5-« 

O 

4-’ 

U 

O 

Pi 



M-; 

O 

o 

LO 



o 

p. 

to 

a 

Di 

P- 

o 

4-> 

CO 



27 



the early system response for all three configurations are 

identical. On the other hand, the configurations of 1 , 4 and 

5 have different W values. As understandable from 

max 

Figure (3.10), the different h' give much effect to the 

in 3. X 

early time response. So, the control-rod characteristic 
becomes the important factor in determining the system 
response for the early values of time. In this paper, the 
control-rod characteristic of 1, which has an application 
in a particular reactor study [Ref. 2] , is used in the 
following sections. In addition, a faster rod response 
characteristic is used in Section D where maximum reactivity 
insertion is investigated. 

C. EFFECT OF TEMPERATURE COEFFICIENT 

The temperature coefficient is assumed to be a constant, 
but in a practical reactor it may vary over some range. 

To investigate the effect of some range of the tempera- 
ture coefficient, the following temperature coefficients 
are used and the responses are shown in Figure (3.11). 

= -8 X 10 ^ 

= -3.5 X lO'^ 

= 1 X 10 ^ 

= 10 X 10'^ 

As the figures indicate, for values of the temperature 

_ 3 

coefficient about the nominal value of -3.5 x 10 , these 

temperature coefficients do not have much influence to 
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Figure (3.11) 

Effect of Temperature Coefficient 



the system response. It works to increase or decrease 
the overshoot. 

D. CONSTRAINT OF REACTIVITY INSERTION 

Total reactivity 6k which is the sum of the reactivity 
due to the control rod, and the temperature reactivity is 
expressed as: 

6k = 6k^ + 6k^ (3.1) 

In the previous sections, the magnitude of the reactivity 
was not constrained. As shoii’n in Section C, the magnitude 
of temperature reactivity is not very influential in the 
system response. The interested point is the constraint 
applied to the control-rod reactivity rate. The movement 
of the control-rod should be restricted to a maximum value 
from physical viev/points. 

Now, the graph of the control-rod reactivity of the 
initial power level 80% is shown in Figure (3.12). Hence, 
as described in Section B in Chapter III, the faster 
control-rod is chosen. Three different constraints of 
control-rod reactivity are considered. These throe con- 
straints are interesting points, because the first works 
to reduce both the first and the second overshoots, and the 
third constraint is entirely below the actual reactivity 
except the starting area. These situations are shown in 
Figure (3.12) . 
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Reactivity 
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Figure (3.12) 

Reactivity of Control Rod for 80 to 100% Power Change 



The Figure (3.13) shows the relation of the system 
response and the control-rod reactivity. The obvious fact 
is that the control-rod reactivity gives much effect, 
especially stability, to the system response. The influence 
of the constraints of the control-rod reactivity is shown 
in Figure (3.14). As shown in Figure (3.14), the 

- 3 

constraint of 6k^ = 3.9 x 10 is the most effective one, 
and the desirable system response of fast rise and minimum 
overshoot is obtained. So the suitable constraint on the 
reactivity insertion of the control rod contributes to the 
stability of the system while improving the system 
performance . 
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Constraint 
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Figure (3.14) 
Influence of Constraints 



IV. OPTIMiXL CONTROL PROBLEM 



There is a possible need to investigate the application 
o£ optimal control theory to nuclear reactor system in 
order to improve the system performance. In this thesis 
the realization of optimal control systems will be investi- 
gated for step demand changes in power. The procedure is, 
however, applied for any deterministic demand input. 



A. OPTIMAL CONTROL LAW 

The optimal control theory has been derived in many 
literatures [Refs. 5, 6 and 7]. The optimal control problem 
is to find a control y*ell which causes the nonlinear system 

x(t) = f(x(t), P(t), t) (4.1) 



to follow a trajectory x*cX that minimizes the performance 
measure 

tf 

J = h(x(t.), tr) ^ f ^ g(x(t), p(t), t) dt (4.2) 

^0 



For most cases of practical viewpoint, such as the control 
of a nuclear reactor, the optimal control must be realized 
in a feedback manner. 



y = f (x(t) , t) 



(4.3) 



while the mathematical solution of this problem is established, 
it involves the solution of a nonlinear partial differential 
equation . 

Optimal -controller synthesis is complicated by the 
nonlinear model which is utilized to approximate reactor 
dynamics. For linear systems, however, it is possible 
to obtain analytical solutions to such problems. Further, 
it is shown in this section that the synthesis of the control 
which minimizes an integrated quadratic in state and control 
only requires a linear feedback of state variables. 

Assume for convenience that the dynamical behavior of 
the system is approximated by the linear model 



The problem is to find the admissible p(t), or, preferably, 
y(x), which controls equation (4.4) with respect to some 
reference x*(t) so as to minimize 



x(t) = A(t) x(t) + B(t) p(t) 



(4.4) 



where x(t^) = x^ 



t 



f 



J = / 
t 



y(t), t)dt 



(4.S) 



o 



v;here 




+ y* ♦ R(t) y] 



(4.6) 
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and R(t) is a real symmetric positive-definite matrix and 
Q(t) is a positive-semidefinite matrix. 

Obviously, this problem may be solved by the maximum 
principle (or Hamilton's equations), or by dynamic program- 
ming (or the Hamilton- Jacobi-Bellman equations). To use 
the Hamilton- Jacobi -Bellman equation, the Hamiltonian is 
formed : 

H(x(t) , u(t) , Iy ’ 



ajT 

+ M . [A(t) x(t) + B(t) y(f)] (4.7) 

A necessary condition for y(t) to minimize H is that 
^ = 0; thus 

If hCt), ^ , t) 

= R(t) u(t) + -11 = 0 C4.8) 



Since the matrix 

2 

|-| = R(t) (4.9) 

y 

is positive definite and H is a quadratic form in y. the 
control that satisfies equation (4.8) does minimize H. 
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Solving equation (4.8) for vi*(t) gives 



q*(t) = 

which when substituted 
H(x(t), y*(t), ^ , t) 



-R'ht) II 



in equation (4.7) yields 
_ 1 T 1 -1 



^ aj , 

3x ~ ~ 



(4.10) 



T ^ 
9x 



(4.11) 



the Hamilton- Jacobi-Bellman eauation is 



9 J(x,t) 

TE 



1 rp 1 9J(X,t) _ -I p 9J(Xjt) 

+ ^ x' Q X - ^ B R ^ b' 



9.x 



9J(x,t) 

"9x 



A X 



(4.12) 



From equation (4.5) the boundary condition is 



J (x(t^) ; t^) = 0 



(4.13) 



Since the minimum cost for the discrete linear regulator 
problem is a quadratic function of the state, it seems 
reasonable to guess as a solution of the form [Ref. 5] 



J (x(t), t) = I J (t) K(t) x(t) (4.14) 

where K(t) is a real symmetric positive-definite matrix 
that is to be determined. Substituting this assumed solution 
in equation (4.12) yields the result 



38 



IT* IT IT -IT T 

•^-x^Kx + ^jX Qx-^^x KBR B^Kx + x KAx=0 (4.15) 

where generally all these matrices and vectors are functions 
of time. Equation (4.15) must be valid for all x(t). 

Using the matrix property, this equation becomes: 

K(t) + K(t) A(t) + A'^(t) K(t) 

-K(t) B(t) R'^(t) B'^(t) K(t) + Q(t) = 0 (4.16) 

and the boundary condition is 

K(t^) = 0 (4.17) 

Once K(t) has been determined, the optimal control law is 
given by 

y*(t) = -R‘^(t) B'^(t) K(t) x(t) (4.18) 

If the system is completely controllable, it is sometimes 
desirable to let t£ so that the optimal control 

accurately maintains the desired terminal state once it is 
reached. Further, Kalman [Ref. 8] has shown that if 

oo 

J = \ f (x'^(t) Q x(t) + p'^(t) R y(t))dt (4.19) 

where Q and R are positive-definite, constant, symmetrical 
matrices and the system [equation (4.4)] is time invariant 



39 



(A and B are constant matrices) , then 

lim dK(t) 
t->«> dt 

In this case, equation (4.16) becomes 

ka + a'^k-kbr‘^b'^k + q= o 

For this problem, the optimal control is 

p*(x) = -R'^ b’^ K X = - X 

and the synthesis requires a constant linear feedback of all 
the state variables as shown in Figure (4.1). 



(4.20) 



(4.21) 



(4.22) 




Controller 

Figure (4.1) 

Optimal Feedback Controller for Linear Regulator 
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B. APPLICATION OF A SECOND -ORDER MODEL 

The second-order differential equation is expressed as 






(4.23) 



This state equation is: 



• N 



2 



-^0 -^1 





f \ 

X-I 




'o ' 




1 


+ 










b. 




2i 




O' 



U 



(4.24) 



The elements of K are determined by solving equation (4.21) 
The corresponding expressions for the optimal feedback 
parameters are obtained by using equation (4.22), and 
developed in Appendix B with the result. 

1 



^0 



0 



'1 u2i2 






(4.25) 



^1 = 






b„ l^^l r ^^0 ^^0^0^ 



(4.26) 



where 





^^1 


0 ' 








'^11 


^2 


= 






> 


? = r ’ 


II 

i 








0 


^ 2 , 








,^12 


^22 



and 



, S'^22 

^0 r ’ ^1 r 
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and selection of the positive parameters corresponds to the 
positive definite requirement on k. 
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V. OPTIMUM APPROXIMiMION OF NUCLEAR 



REACTOR SYSTEM BY A SECOND- ORDER MODEL 

The nonlinear higli-order system described in Chapter II 
which is written below ' 



x=Ax+Be (5.1) 

where ’ 

e = X* - X 

and y^^^=c"^x (5.2) 

may sometimes be approximated with sufficient accuracy 
by a linear model for small variations in state and control 
To apply the optimal control law to the nuclear reactor 
system, this optimum approximation problem is the most 
necessary thing. This idea has been developed by Sinha and 
Bereznai [Ref. 9]. But in this paper, the different 
error criteria are used. 

A. DESCRIPTION OF METHOD 

Considering a discrete set of values of Yr.>. taken over 

~ v.t j 

a suitable interval of time. 
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Y = {y^, 72, . . . , y^> (5.33 

where y. = 

This set represents discrete sample points of the response 
of the system, which is obtained by solving equation (5.13, 
using the digital computer. These results are shown in 
Chapter III. The continuous output response y^-^^ is sampled 
at sufficiently close intervals of time lest the significant 
information is lost. The objective is to find another 
output set y* in optimal fashion, associated with a second- 
order model described by the equations. 

X* = A* X* + B* p (5.43 

y=^ = C X* (5.53 

such that, for the same input, the following objective is 
satisfied . 

A scalar error performance function J is minimized 

J = g['d (>’i - ypl C5.6) 

which is some .suitable function of the errors (y. - y?3 with 
a vector weighting factor attached at each sampling 
ins tant . 

The choice of the error criterion expressed by 
equation (5.63, has a direct effect on the parameters of 
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the approximating model. Since the purpose o£ the error 
criterion is to measure the extent to which the response of 
the model deviates from that of the actual system, the main 
problem is how to express this deviation numerically. It 
is normal practice to consider the norm of the output error 
and raise it to some power p. This objective becomes to 
minimize a summation of deviation. 

J = E W. I J y. - y* I 

for the single-output case, v.'here is a scalar sequence, 

and the weighting factor is unity, equation (5.7) means 
a measure of the area between the curves when P=l, and the 
mean square error when P=2. 

Response 




and the Optimal Second-Order Model 



45 



Consider the system response the model response 

^*(t) shown in Figure (5.1). It is required to find a 

model of the system so that the error criterion is mini- 
mized. The following two error criteria as the performance 
measure are used to approximate the response yj-^^ to the 
response of some model. 

(1) minimization of 



I CXi 
1 



yp 



(5.8) 



(2) minimization of 



= E {(y. - v*)^ + (^^ 
^ '1'^ ^dt 



£ 1 ) 2 , 

dt ^ 



(5.9) 



by considering the derivative of y^ and y|, the more signifi 
cant measure is attained v;here the slope is changing rapidly 

B. OPTIMUM APPROXIMATION USING PATTERN SEARCH 

The problem of approximating high order system by a 
second-order model in an optimum manner can only be solved 
for specific error criteria by the presently available 
techniques . 

Since an analytical solution to this problem does not 
appear feasible, various search techniques may be considered 
While gradient methods are generally quite efficient in 
finding a minimum, the necessity of finding the partial 
derivatives of arbitrary error criteria with respect to all 
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the model parameters, becomes a great disadvantage. 

On the other hand, direct search techniques involve 
evaluating the effect of sequential parameter changes in 
an organized manner. The pattern-search technique of Hooke 
and Jeeves [Refs. 10 and 11] was selected as a suitable 
method for this task. Appendix C shows an outline of the 
pattern- search technique, which is contained within the 
IBM 360/67 library under the subroutine name of DIRECT. 

C. INITIAL GUESS OF COEFFICIENTS OF A SECOND-ORDER MODEL 
An important problem associated with search technique 
is the selection of the starting parameters, since these 
have considerable influence on the convergence of the 
process, and on the probability of locating a local 
minimum of a performance measure. Therefore, the reasonable 
initial guess to compute the parameters using subroutine 
DIRECT is needed. The characteristics of the second-order 
linear system have been developed in many manners [Ref. 12] . 

Since the response of a nuclear reactor system is approxi- 
mated by a second-order model, the characteristic of that 
is applicable to determine the coefficients. The differen- 
tial equation of a second-order system in the form of the 
natural frequency, W^, and the damping ratio, c, is expressed 
as 

2 

^ + 2 C h ^ X = KU(t) (5.10) 

j. 2 ^ndtn 
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The response x(t) and its first derivative arc obtained by 
standard methods, and are shown below. 



x(t) = 



AK 

W ' 
n 



[1 



/ 



- cV,' t . 
e n sin 



(W *^1 - 
^ n 






(5,11) 



1 - C2 



' r.-. AK e 

x(t) = ^ sin (W^ 1 - C t) (5,12) 

n -'i _ 

Figure (5,2) is a plot of x(t) and x(t). Since the input 
is a unit step, the value of A equals to one. To make the 
steady state value of x(t) unity. 



x(t ) 





(5.15) 



2 

Therefore, if K = , the steady state value of x(t) 

becomes unity. 

The value of x^^^^ may be found by a setting equation 
(5.12) to zero, finding t^, and substituting this value of 
tp into equation (5.11) 



' /-.u ^ n AK 
^ n / 



-CW t / =• 

^ sin (W t) 

2 

1 - C 



(5,14) 



Then 



W t 
n p 




(5,15) 
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/ 




and 



Figure (5.2) 
x(t) and x(t) vs. time 



x(t ) 
P 



X 

max 



CTT 




(5.16) 



for a special case of a second-order system expressed in 
equation (4 . 23) , 



^i ^2 



^2 " ^0 ^1 " ^1 ^2 



O 

As stated previously, choosing a^ = = 2?W^ 

and U(t) = unit step input, the steady state value of x^ is 
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x^(t ) 



(5.18) 




The characteristic response of the reactor system of the 
initial pov^;er level 90% solved by digital computer is shown 
in Figure (5.3) 




Figure (5.3) 

Response of Nuclear Reactor System (90-100%) 

Substituting t =2.48 and x = 1.0144 into equations 
® p max 

(5.15) and (5.16), the computed natural frequency and 
damping ratio are; 

C = 0.8035 

W =2.128 
n 

Then, a. = W ^ = 4.528 

’ 0 n 

a, = 2cW = 3.42 
1 n 



50 



Using these values as the initial guess coefficients, the 
computed optimal parameters by search routine are: 

ap = 1.554 

a^ = 1.51 

There is considerable difference between the initial guess 
values and the computed optimal coefficients. In a practi- 
cal point and considering the nonlinear nature of the system, 
these initial guess values are sufficiently near the 
optimal coefficients. Since there is a definite necessity 
to compute initial guess values, this method is shown to 
be acceptable . 

D. COMPUTER PROGRAM OF A PATTERN SEARCH 

A computer program has been v;ritten that uses a pattern 
search subroutine to find optimum second-order models for 
high-order systems. The program has the following summarized 
features . 

(1) The output responses of the nuclear reactor system to 
unit step function are available at discrete uniform 
intervals of time. 

(2) The parameters of the model to be used for starting 
values are given by the user. 

(3) The program assumes a uniform weighting sequence. 

(4) The second-order model is solved by the Runge-Kutta- 
Gill fourth-order method. All calculations are in 
double precision. 

(5) The objective functions to be minimized are 
equations (5.8) and (5.9). 
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E. MANIPULATION OF THE RESULT 



Using the search routine, the coefficients of the 
approximated second-order model are manipulated with the 
results shown in Table (5.1), and the plot of that is 
Figure (5.4) . 

The difference in the model coefficients using the 
performance measure and that of the performance measure 
is small. This is clearly shown in Figure (5.5). 

The system responses and the approximated second-order 
model of the performance measure J 2 to various initial po^^^er 
levels are shown in Figure (5.6) through Figure (5.11). 

F. CONSIDERATION OF A THIRD-ORDER SYSTEM 

Comparing the system response with the second-order 
model as shown in the previous part, there is some appre- 
ciable difference. But, as tiie considered system has a 
ninth-order nonlinearity, to approximate that by a second- 
order model, the difficulty exists in the natural sense. 

So this difference is considered acceptable. As a trial, 
even though the optimal control law of the third-order model 
is not proposed, the approximated third-order model is com- 
puted using the same search routine. The equation of that 
is : 

?? ^ d? ^ dT " ^3 (5.19) 
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Table (5.1) 

Coefficients of Approximated Second-Order Model 
to Various Initial Power Level 



Initial Power Level 

Nic = 0.95 (28 sample points) 

= 8.191 X 10'^ 
a^ = 4.399 
a^ = 2.149 

= 1.399 X lO'-^ 
a^ = 4.008 
a^ = 1.809 



N 



IC 



0.9 



J 



1 



J 



2 



(28 sample points) 
6.60 X 10'^ 
aQ = 1.651 
a^ = 1.614 

5.264 X 10'^ 
a^ = 1.554 
a^ = 1.509 



Nic = 0.8 (28 sample points) 

= 3.792 X 10'^ 
ap = 0.507 

a^ = 0.950 

= 1.539 X 10‘^ 
aQ = 0.478 
a^ = 0.885 



Time Interval 
0.25 



0.25 



0.25 
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r 




Njc = 0.6 



J, = 1.791 X 10 



(28 sample points) 
-2 



0 . 5 



= 0.14 
= 0.511 



= 3.529 X 10 



-2 



a^ = 0.137 
a^ = 0.496 



Nic “ ^ 



J, = 3.705 X 10 



(28 sample points) 
2 



0.5 



Oq = 0.0804 



a^ = 0.358 



J., = 6.463 X 10 



-2 



aQ = 0.0795 
a.^ = 0.351 



Njc 



(28 sample points) 



J. = 



Jo = 



7.424 X 10’^ 




ag = 2.135 X 


lO’^ 


a^ = 5.689 X 


lO'^ 


1.409 X 10'^ 


- 2 


ag = 2.282 X 


10 


a^ = 7.353 X 


10’^ 



0.5 
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0.2 



0.5 



N 



IC 



(28 sample points) 
= 6.097 X lO'^ 

Eq = 1.564 X 10'^ 

= 6.097 X 10'^ 

= 3.560 X 10‘^ 

Eq = 1.633 X 10'^ 

a, = 6.966 X 10 ^ 



5 5 




Power Level 



Figure (5.4) 

Coefficients Plot of Second Order Model 
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Figure (5.6) 

Comparison of Step Response of System 
and Model to 95% Initial Power Level 




Comparison of Step Response of System 
and Model to Initial Power Level 80% 



/ 
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Comparison of Step Response of System 
and Model to Initial Power Level 60% 




Comparison of Step Response of System 
and Model to Initial Power Level 50% 
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Comparison of Step Response of System 
and Model to Initial Power Level 20 % 
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where U(t) = unit step input. The obtained data of the 
optimal third-order system to initial power level 50% is: 



= 6.213 X 10 
a^ = 2.336 
a^ = 8.10 
= 10.5 
a.. = 10.64 
= 7.7 X lO'^ 

aQ = 2.16 

a^ = 8.53 
3 . 2 - 10.63 
a^ = 10.88 

The clarified plot is shov.Ti in Figure (5.12). Also, the 
second-order model is added in this graph. It is observed 
that the third-order model is m.ore realistic than the 
second-order model to approximate the actual system. Since 
it is desirable to construct the optimal control law of 
the third-order model, that work will be expected in the 
future . 
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3rd Order 




c\o 
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Figure (5.12) 

Approximated Third-Order Model to System Response of Initial Power Level 50 



VI. SUBOPTIMAL CONTROL USING 



A SECOND ORDER MODEL 



The suboptimal control of a reactor system using optimal 
second-order model is made possible by the fact that the 
state variables available for feedback are the actual output 
power level and its rate of change. The optimum linear 
feedback of these variables for the case of the second- 
order model with an integral quadratic cost function exists, 
and the appropriate feedback coefficients can be determined 
explicitly in terms of the parameters a^, a^^ and the 
constants used in the cost function. These realize a sub- 
optimal control for the reactor system. The block diagram 
representation of the suboptimal controller, the reactor 
system and the optimal controller for the second-order model 
are shown in Figures (6.1) and (6.2). 

Using Table (5.1), the calculated k^ and k^ are shown 
in Table (6.1). Hence, the weighting factor = q 2 = r in 
equations (4.25) and (4.26) is used to compute kQ and k^ . 

In this system and model with the suboptimal controller 
and the optimal controller, one modification is needed 
to make the steady state value unity. 

From Figure (6.2) 




a 



0 



( 6 . 1 ) 



TTPT 
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Figure (0.1) 

Reactor and Suboptimal Controller for Step Demand Changes 
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Figure (6.2) 

Second-Order Model and Optimal Controller for Step Demand Chang 



Table (6.1) 
Controller Parameters 
to Various Initial Power Level 

Initial Power Level 



95% 


0.4142 


0.736 


90% 


0.4142 


0.602 


oo 

o 


0.4142 


0.631 


60% 


0.4142 


0.869 


50% 


0.4142 


1.144 



and 

U(s) = k U^(s) - (kp + k^s) x^(s) (6.2) 



Substituting (6.2) into (6.1) 



x^(s) a^k 

+ (a^ + k^) s + a^ (1 + k^) 



Using the final value theorem 

V = lira 

^l(oo) s->0 



s x^(s) 



k 









(6.3) 



(6.4) 



Therefore, the feed forward controller of the amplitude 
k = 1 + kg to make the steady-state value unity is needed. 
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To make the comparison possible, between the responses 
of the suboptimal controller and the corresponding optimal 
model controller, and that of the system and the model 
without any controller. Figure (6.3) shows a plot of the 
system response to initial power level 90% without the 
controller and wit^i the controller. The describable fact 
is that the maximum overshoot of the system response with 
the suboptimal controller is reduced. Also, Figure (6.4) 
shows a plot of the model response to the same initial power 
level without the controller and with the controller. The 
desirable reduction of the time to reach the 100% power 
level is observed in the model response with the optimal 
controller. So, the considered suboptimal controller 
works to make the reduction of the maximum overshoot of the 
system response with no effect on the early time response. 

On the other hand, the optimal controller of the model 
does affect the early time response by making it faster. 

Table (6.2) is the comparison of the overshoot and the time 
to reach the first full power (100%) of the system and the 
model, without controller and with controller. 

The responses between the system and the model to 
various initial power level with the controllers as shown 
in Figure (6.5) through Figure (6.9). Comparing these 
with Figure (5.5) through Figure (5.9), the controllers work 
to effect the model rise time in such a manner to make the 
overshoot of the system and the model occur at the same. time, 
and have approximately the same value. 
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System 103.81 102.30 59.6% 7.99 8.02 +0.03 

Model 102.42 105.14 10.23 7.59 




Figure (6.3) 

Comparison o£ System Response of Initial 
Power Level 90% without Controller 
and with Controller 




Figure (6.4) 

Comparison of Model Response of Initial 
Power Level 90% without Controller 
and with Controller 
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Power Level (I) Power Level (%) 



Model 




Figure (6.5) 

Step Responses of System and Model 
with Controllers to Initial Power Level 95% 




Figure (6.6) 

Step Response of System and Model 
with Controllers to Initial Power Level 90% 
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Step Responses of System and Model 
with Controllers to Initial Power Level 80 




Step Responses of System and Model 
with Controllers to Initial Power Level 60 
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Initial Power Level 50% 



VII. REALIZATION OF ADAPTIVE CONTROLLER 



The effect of the suboptimal reactor controller and the 
optimal model controller is described in the previous 
chapter. There the controller parameters are fixed. It 
is possible to design programmed time variations of 
controller parameters to achieve instantaneous near optimum 
control for the system at all times. In addition, the 
system parameter such as the temperature coefficient 
change \\rith time and the existence of such conditions is 
the reason for application of an adaptive control system. 
Eveleigh [Ref. 13] has defined the adaptive control as 
follows : 

"If an index of performance (IP) is available which 
indicates the system's instantaneous or short-term 
average performance quality, and if a control loop 
is set up to optimize the IP automatically by adjusting 
controller parameters, the parameter-adjustment configura- 
tion is called an adaptive control loop. It is important 
to note that an adaptive controller is a parameter- 
adjustment loop above and beyond the normal feedback 
used to control position, velocity, and the like. 

Adaptive control is thus an effort to extend basic 
optimum-control concepts to time-varying systems." 

A. PROCEDURE OF ADAPTIVE CONTROL 

Hence, the concept of the adaptive control, which is 

shoiv^n in Figure (7.1) as used in Reference 2, is employed. 

This procedure is described below. 

The model parameters are precomputed, over the expected 

range of demanded power level changes, to be used as starting 

values for the optimization. As the model coefficients 
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Figure (7.1) 

Block Diagram o£ Adaptive Control System 
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4 




are updated, the parameters of the feedback controller are 
computed from equations (4.25) and (4.26), depending on the 
nature of the cost function. Here, also, the weighting 
factor = q 2 = r is used to make the comparison possible 
in the previous chapter. That is, the cost function has the 
following form: 

CO 

J = y / (x^ + U^) dt (7.1) 

z 0 ~ 

Initially, the nuclear reactor is assumed to be operating 
at a steady state power level. At time t = 0, a step 
change to a new operating level is required. Using the 
model parameters appropriate for the demanded change and 
the desired cost function, the controller parameters are 
evaluated, and the system begins the transition to the new 
power level. Over several sampling intervals, the value 
of the system response at each sample point is stored. 

This means neither updating of the model nor recomputation 
of the controller parameters takes place. 

The fitness between system and model depends on the 
accuracy of the starting model. The model is reoptimized on 
the basis of the observed system response. The next step 
is to recompute the controller parameters. After completing 
these procedures, the most optimized controller parameters 
of some interval is found. This optimized controller 
introduces the desirable system response. Repeating these 
steps, the system response of entire intervals is obtained 
in the most optimal manner. 
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As described in Section B of Chapter III, the high 
gain of the control rod rate gives the response time 
much influence. To make faster response, the following 
control rod is used instead of equation (2.5). 




Figure (7.2) 

High Gain Control Rod Characteristic 

The system response without controller of initial 
power level 50% using the high gain control rod is shown 
in Figure (7.3). Also, the approximated second order 
model is shown in the same plot. The time to reach 100 
value of the system is earlier with the higher gain rod. 

Apparently, judging from this graph, the second-order 
model does not fit to the system response. Though the 
third-order model describes the system more accurately 
than the second-order model as shown in Figure (5.12), 
the usage of the third-order model is not established 
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Figure (7.3) 

Responses of System of High-Gain Control Rod 
and Approximated Model without Controller 



theoretically, and therefore is not applied to this 
problem. 

Figure (7.4) shows a plot of the system response with 
the suboptimal controller to initial power level 50^ and 
that without controller to the same initial power level. 
The suboptimal controller based on the optimal control 
law of the second-order model has the great effect to 
reduce the overshoot and the system response reaches 
steady-state value quickly. Apparently, the usage of the 
linear feedback controller is recommended to control the 
nuclear reactor to reduce the large overshoot and to 
reach the steady state faster. Figure (7.5) is a graph 
of the responses of the system and the model with the 
controllers. Comparing this with Figure (7.3), the 
second-order model with optimal controller more closely 
represents the reactor system than when no controller is 
used . 

B. APPLICATION OF ADAPTIVE CONTROL 

As described in the previous section, the procedures 
of the adaptive controller are applied by using an obser- 
vation interval of 15 seconds duration and the adaption 
intervals of 0.25 second. The obtained parameters 
that are computed repeatedly until sufficient accuracy is 
attained, are shown in Table (7.1). 

The corresponding responses of system and model with 
the adaptive controller are shown in Figure (7.6). Also, 
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System without Controller 
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Figure (7.4) sec 

Effect of Suboptimal Controller 
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Figure (7.5) 

Step Responses o£ System and Model with Controller 
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Figure (7.6) 

Step Responses of System and Model \\rith Adaptive Controller 



Parameter Changes during Power Level Transitions 
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Figure (7.7) shows the system responses with the non- 
adaptive controller and with the adaptive controller. A 
comparison of Figure (7.5) and Figure (7.6) shows that 
there is some improvement between the model response 
with the optimal controller and that with the adaptive 
controller. In Figure (7.7), there is no appreciable 
difference between the system response with the non- 
adaptive controller and that with the adaptive controller. 
This means the concept of adaptive controller is not 
effective as long as the system parameters are constant 
and unity weighting factors are used in the cost function. 

As a reference, the response of the second-order model 
used to develop the adaptive controller is shown in 
Figure (7.8) with the system response using the adaptive 
controller. 

To show the effect of adaptive controller 
clearly, another example is studied. The theory of 
adaptive control is applied to the system of initial pov^er 
level 801 shown in Figure (7.9). The tabulated result 
is shown in Table (7.2). 

Figure (7.9) is a plot of the system responses with 
the non-adaptive controller, and with the adaptive con- 
troller. Though there is very little difference between 
them, they are almost considered identical. 
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Figure (7.7) 

Comparison of System Responses with Adaptive 
Controller and Non-adaptive Controller 
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Figure (7.8) sec 

Responses of System and Updated Second-Order Model 
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Figure (7.9) 

System Responses with Suboptimal Controller 
and with the Adaptive Controller 
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VIII. INFLUENCE OF COST FUNCTION 



As having been described, the cost function is: 



00 

J = i / (x^ Q X + rp^) dt 

z 0 ~ 



( 8 . 1 ) 



In the previous chapter, the controller parameters were 
computed for the case that Q is an identity matrix and r 
equals to unity. 

In this chapter, the various weighting factors for 
the cost functions are considered, depending on the optimal 
control law. The selection of appropriate values for Q 
and r is important, and goes to the core of the optimal 
regulator problem. The process of selection consists of 
assigning certain values to Q and r, incorporating the 
resulting controller into the system. 

Equations [4.25 and (4.26) of the controller parameters 
are rewritten for convenience 



^0 



^1 = 






4 ro2 4 . ^1 k2i2 

t— [a^ + — b.J 



^ 1.1 r2 






where q^ and bhe components of a matrix q, r is the 

relative weighting factor. 

These weighting factors of the cost function decide 
the controller parameters. 
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The following values are tried to check the behavior 



of the actual system. 




0.1 



The computed controller parameters are: 



kp = 0.0488 
= 0.1426 




kQ = 0.4142 
k^ = 0.797 




k^ = 2.317 
k^ = 3.008 




k^ = 3.583 
k^ = 4.558 



These responses are shown in Figure (8.1). The stability 
of the system with the suboptimal controller depends on 
weighting factors of the cost function apparently. Since 
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Figure (8.1) 

System Response as the Function of the Cost Function 



the optimal control law is based on the quadratic cost 
function, it is essential to select the most reasonable 



cost function, 
resulting fact 



as described in Reference 14. The 
is consistent with Reference 14. 
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IX. CONCLUSIONS AND RECOMMENDATIONS 



The following conclusions and recommendations are 
offered as a result of the investigation contained in 
this report. 

A. CONCLUSIONS 

1. The simulation of the dynamic behavior of the 
nuclear reactor system on the digital computer to make 
feasible this investigation was obtained lAfith some 
difficulty. 

To obtain the stabilized system response required 
considerable effort, because the system is nonlinear and 
described by a ninth-order differential equation. The 
unknown factors such as the specified values of the 
nuclear reactor parameters complicated the difficulty of 
implementation of the system equations. 

2. In this paper, the gain of the control rod, the 
reactivity-rate constraint and the effect of the tempera- 
ture coefficient have been studied. The high gain of the 
control rod causes the faster response, and the reactivity- 
rate constraint suitably chosen makes possible the more 
desirable response regarding the stability, particularly 

to reduce the deviation from the desired power level. The 
magnitude of the temperature coefficient has a lesser 
effect on the system response. 
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3. It has been shown that the results from the 
optimal second-order model of the reactor may be used 
effectively for suboptimal control of the reactor system. 
This suboptimal controller of the reactor reduces the 
overshoots of the system response. Considering these 
overshoots are not desirable practically, its usage 
provides an improvement in system performance. To find 
the optimal second-order model, the computer program of 
search technique is applied to obtain the second-order 
model parameters. 

4. An adaptive control system, which may be realized 
in practice to control the power level changes of a 
nuclear reactor, has been proposed. 

The parameters of the second-order linear model 
are continuously updated so that the model accurately 
represents the behavior of the system. Also, the parameters 
of the controller are updated to achieve the concept of 
the adaptive control. The result of this study shows no 
improvement of the system response using the adaptive 
controller for the operating condition considered. The 
adaptive control, therefore, would not have to be on- 
line, and could be used to update the model parameters 
to account for slowly changing parameters of the reactor 
system. 

5. The weighting factor of the cost function gives 
the stability of the system response much influence. 
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Based on physical considerations, the weighting factor 
should be chosen as correctly as possible to provide the 
desired response and this would be a trial and error 
procedure . 

B. RECOMMENDATIONS 

Further workers in this area should consider the 
following facets of nuclear reactor systems. 

1. The nuclear reactor system studied in this report 
did not consider the noise problem. For the purpose of 
the present paper, it was assumed that an instantaneous 
noise-free measure of the reactor power level is available. 
However, as the temperature channel signal indicates the 
power level of the reactor at an earlier instant, because 
of the finite transport time between the reactor core 

and the temperature transducer, this signal has a con- 
siderable noise component due to the turbulent coolant 
flow. 

Therefore, the noise consideration is strongly 
recommended for future studies. 

2. Since the proposed scheme relies on the use of 
search routines, improvement in the efficiency of the 
method particularly as applied to digital process computers, 
would greatly enhance the usefulness of the new technique. 

3. In this report, the effectiveness of the adaptive 
control has not been shown. Consideration should be 
given to the adaptive controller for different weighting 
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factors in the cost function to find out whether an on-line 
adaptive scheme would be required. 

4. In much of optimal control, the natural formulation 
is by state variables. A w'ell-known and popular result is 
that for the linear regulator problem with quadratic 
performance index as applied in this study. It has been 
shown [Ref. 15], however, that the solution is far from 
realistic or optimum in an engineering sense. Nevertheless, 
the results and scheme obtained for the linear quadratic- 
index regulator problem have been forced indiscriminately 
on the linear time- invariant problem. 

Since the weakness of the s tate -variable formula- 
tion in coping with linear feedback control system design 
is especially apparent in the sensitivity problem, which 
is one of the primary reasons for the use of feedback, the 
sensitivity theory by Horowitz and Shaked [Ref. 16] should 
be investigated. 
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APPENDIX A 



NUNERICAL VALUES OF THE REACTOR MODEL PARAMETERS 
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APPENDIX B 



DERIVATION OF EQUATIONS ON AND 
The equation is : 



ka + a'^k-kbr'^b^k + q = o 



where 



A = 


f 

0 


1 


B = 


/ \ 
0 




, "^0 






. ^0, 



9 = 


'^1 


0 ' 


R = r , K = 


’^11 


^12 




0 






^^21 


^22 



Computing term by term 
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0 0 
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-IT 0 

K B R K = — - 



^12^21 



^ 22^21 



^ 12^22 



'22 



After substituting these matrices into the equation, 
rearranging, the following four equations are obtained 

u 2 



-a„ k, ^ - a„ k 



0 12 0 12 r '"12 ^1 



° k^- + q = 0 



a^ k^^ + k, , - a, k 



0 



k,. k.. = 0 



0 ^"22 11 1 ^"12 r ^"12 ^"22 



k,, - a, k^- - a^, k 



k,o k,, = 0 



11 1 21 0 ^"22 r *"12 ^"22 
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12 1 22 ^"21 1 22 r 22 ^2 



° + q = 0 



From the first one of the above four equations 
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From the optimal control law, 



U* = - R'^ K X = -k^ X 
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APPENDIX C 



OUTLINE OF THE PATTERN SEARCH 

As devised by Hooke and Jeeves [Ref. 10] , pattern 
search is a direct search method which results in relatively 
efficient search along straight ridges or ravines. It is 
attempted with pattern search to establish the pattern 
of successful search points in the immediate past from 
which plausible future search points are predicted. The 
method for the two-dimensional case is illustrated in 
Figure (C.l]. The search progresses from an initial base 
point x°, and a minimum of f(x) is to be found without the 
use of derivatives of f(x). A small displacement d^^ 
from x° in the x^ direction is effected, and f(x^ + d^, xp 
is evaluated and compared with f(x|, • If the latter 

is smaller than the former, a small displacement in the 
-x^ direction is effected, and f (x| - d^ , xp is compared 
with f (xp ^p • If the latter is greater than the former, 
a small displacement d -2 is made from (xj - d^, xp in the 
X 2 direction, and f(x^ - d^ , x°) is compared with f(xj - d^ , 
^2 * ^2^ ‘ Assuming the latter is the smaller of the two, 
a pattern move is made in the direction established by the 
preceding successful moves to the point x^ in Figure (3.1]. 
At x^ , the small moves which proved to be successful around 
x° are repeated and if successful lead to a second pattern 
move to X . From x , small displacements of magnitude d^ 
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Figure (Cl) 
Pattern Search 
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along the coordinate prove to be unsuccessful, whereas 

2 2 2 2 

the displacement 62 yields £(x^ , x^ + d2) < f(x^ , ^2) . 

A pattern move is made therefore along the line which 
passes through the points (xj^ - d^ x^ + d2) and 
(x^ , X 2 + d2) . 

Note that the step size of the pattern search is made 

progressively larger as the general directions dictated 

by the pattern prove successful. When failure occurs, 

£ 

as at point x , the step size of the pattern move is 

3 3 

reduced as depicted, to point x . At x , the displacement 
3 

d^ from X along the x^ coordinate fail to improve f (x) , as 
do displacement d2 from x along the x coordinate. Reduc- 
tions in d^ and d2 are therefore required at this point, 
and the process starts anew. Ultimately both the dis- 
placement and the pattern step sizes are reduced below 
some preassigned limit, and the search is terminated. 
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